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Abstract 

This article studies the viscous flow and heat transfer over a plane horizontal surface stretched non-linearly in two lateral 
directions. Appropriate wall conditions characterizing the non-linear variation in the velocity and temperature of the sheet 
are employed for the first time. A new set of similarity variables is introduced to reduce the boundary layer equations into 
self-similar forms. The velocity and temperature distributions are determined by two methods, namely (i) optimal homotopy 
analysis method (OHAIVl) and (ii) fourth-fifth-order Runge-Kutta integration based shooting technique. The analytic and 
numerical solutions are compared and these are found in excellent agreement. Influences of embedded parameters on 
momentum and thermal boundary layers are sketched and discussed. 
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Introduction 

The fundamental problem of two-dimensional flow due to 
stretching plane surface, initially discussed by Crane [1], is 
involved in various industrial processes such as metal and polymer 
extrusion, drawing of plastic films, paper production etc. Owing to 
such applications, the researchers have discussed this problem 
under various aspects including suction or injection, variable 
surface temperature, convective boundary condition, mass trans- 
fer, mixed convection etc. The three-dimensional flow due to 
plane bi-directional linearly stretching sheet was first discussed by 
Wang [2] . He found an exact similarity solution of the classical 
Navier-Stokes equations. Later, Lakshmisha et al. [3] numerically 
examined the unsteady three-dimensional flow with heat and mass 
transfer over an unsteady stretching sheet. In contrast to this 
problem, Takhar et al. [4] investigated the three-dimensional flow 
of an electrically conducting fluid due to an impulsive motion of 
the stretching sheet. Ariel [4] derived approximate analytic and 
numeric solutions for steady three-dimensional flow over a 
stretching sheet. Xu et al. [5] provided uniformly valid series 
solutions for three-dimensional unsteady flow caused by the 
impulsively stretching sheet. Liu and Andersson [6] considered the 
heat transfer in three-dimensional flow due to non-isothermal 
stretching sheet. The unsteady three-dimensional flow of elastico- 
viscous fluid and mass transfer due to unsteady stretching sheet 
with constant wall concentration was studied by Hayat et al. [7]. 
In another paper, Hayat et al. [8] described the three-dimensional 
flow of Jeffrey fluid due to stretching sheet. Liu et al. [9] 
firsdy discussed the three-dimensional flow due to exponentially 



stretching sheet numerically. Steady flow of nanofluid past a 
linearly bi-directional stretching sheet through Buongiorno's 
model was examined by Junaid et al. [10]. Sheikholeslami and 
Ganji [1 1] discussed the flow and heat transfer of nanofluid 
between parallel sheets in the presence of Brownian motion and 
thermophoresis effects. 

The literature cited above deals only with the case of either 
linearly or exponentially driven velocity of the sheet. Vajravelu 
[12] numerically discussed the viscous flow due to stretching sheet 
when the velocity of the sheet was assumed to obey the power-law 
distribution, i.e He computed numerical solutions for 

various values of power-law index n. CorteU [13] extended this 
problem by considering viscous dissipation effects and variable 
surface temperature. The steady boundary layer flow of micro- 
polar fluid over non-linearly stretching sheet was discussed by 
Bhargava et al. [14]. Radiation and viscous dissipation effects on 
the boundary layer flow above a non-linearly stretching sheet were 
explored by CorteU [15]. Homotopy analytic solutions for mixed 
convection flow of micropolar fluid past a non-linearly stretching 
vertical sheet were obtained by Hayat et al. [16]. KechU and 
Hashim [17] derived analytic solutions for MHD flow over a non- 
linearly stretching sheet by Adomian decomposition method. 
Hayat et al. [18] used modified decomposition method for tiie 
series solutions of MHD flow of an electrically conducting fluid 
over a non-linearly stretching surface. The impact of chemical 
reaction on the flow over a non-linearly stretching sheet embedded 
in a porous medium was investigated by Ziabakhsh et al. [19]. 
Rana and Bhargava [20] computed numerical solutions for two- 
dimensional flow of nanofluid due to non-linearly stretching sheet 
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by finite element method. Shahzad et al. [21] obtained closed 
form exact solutions for axisymmetric flow and heat transfer when 
the velocity of the stretching sheet was proportional to r^. Partial 
slip (effects on the boundary layer flow past a non-linearly 
permeable stretching surface have been addressed by Mukhopad- 
hyay [22]. In another paper, Mukhopadhyay [23] analyzed the 
flow and heat transfer of Casson fluid due to non-linearly 
stretching sheet. Rashidi et al. [24] derived homotopy based 
analytic solutions for flow over a non-isothermal stretching plate 
with transpiration. 

To our knowledge, the three-dimensional flow due to non- 
Unearly stretching sheet has not been yet reported. It is obvious 
that three-dimensional flows can be appropriate in giving more 
clear physical insights of the real world problem when compared 
with the two-dimensional flows. The present work is therefore 
undertaken to fill such a void. The study also assumes that the 
temperature across the sheet is non-linearly distributed. Introduc- 
ing a new set of similarity variables the boundary layer equations 
are first reduced into self-similar forms and then solved both 
analytically and numerically. It is pertinent to mention that 
computation of either approximate analytic or numerical solutions 
of the boundary layer equations governing the flow and heat 
transfer is often challenging (see [25-33] for details). Attention is 
focused on the physical interpretation of parameters including the 
power-law index n. 

Mathematical Modeling 

Let us consider the three-dimensionalincompressible flow over a 
plane clastic sheet located at z = 0 as shown in the Fig. 1 . The flow 
is induced due to stretching of the sheet in two lateral directions. 
Let UK{x,y)=a{x+yy and Vt„{x,y)=b(x+yy be the velocities 
of the sheet along the x— and y— directions respectively with 
a,h,n>0 are constants (see Table 1). T„(x,y) = Too +A{x+y)" is 
the variable surface temperature where ^ > 0 is constant and Too 
is the ambient fluid temperature. Under the usual boundary layer 
assumptions, the equations governing the three-dimensional flow 



and heat transfer in the absence of viscous dissipation and internal 
heat generation/absorption can be expressed as (see Liu et al. [9]) 

du dv dw 



du du du 8^u 
ox dy dz dz^ 



dv dv dv d^v 
dx dy dz dz^ ' 



dT dT dT d^T 

where u,v and w are the velocity components along the x,y and 
z— directions respectively, v is the kinematic viscosity, T is the 
fluid temperature and a is the thermal diffusivity (see Table 1). 
The boundary conditions are imposed as under: 

u = U„ = a(x+y)",v= Vyf = b{x+y)",w = 0, 
T=T^=T^+A{x+y)"a.tz=0, (5) 
u = 0,v = 0,T-fTao asz->oo. 



f f f f 





I I I I I 

= a{x + yf 




Figure 1. Physical configuration and coordinate system. 

doi:1 0.1 371/journal.pone.01 07287.g001 
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Table 1. List of symbols. 





(xj'.r) Cartesian coordinate systenn 




k thermal conductivity 


M,v,u' velocity components along the a — ,1— ,r — 


directions 


h non-zero auxiliary parameter 


f/„,K|, velocity of the stretching sheet along a — 


and y— direction 


' 1^*^ order derivative with respect to i] 


T fluid temperature 




" 2"^ order derivative with respect to y\ 


r,t, wall temperature 




"' 3^^ order derivative with respect to \] 


ambient fluid temperature 




Greek symbols 


a,h,A positive constants 




V kinematic viscosity 


u Power-law index 




K thermal diffusivity 


f\g dimensionless stream function 




Q dimensionless temperature 


Pr Prandtl number 




}] similarity variable 


Cj\-,Cfy skin friction coefficient along a — and v — 


direction 


X ratio of the stretching rates 


A^M local Nusselt number 




'^zx-.'^zy wall shear stress along a — and v— direction 


^vt, wall heat flux 




p density of the fluid 


Re^,Re^ local Reynolds number along a — and y 


— direction 


fi dynamic viscosity 



doi:l 0.1 371/journal.pone.0107287.t001 

We introduce the new .similarity transformations as follows: 
u = a{x+y)"f'(nX v = a(x+y)"^{nX w 



{x+y) 2 z. 



where Pr = v/a is the Prandtl number and A = b/a is the ratio of 
stretching rate along the y — direction to the stretching rate along 
the X— direction (see Table 1). The above equations reduce to the 
case of two-dimensional flow when 1 = 0. Moreover, when 1=1, 
the equations governing the axisymmetric flow due to non-linearly 
stretching sheet are recovered. When Pr = 1 the solution of /' is 
also a solution of 6. The quantities of practical interest are the skin 
friction coefficients and the local Nusselt number which are 
defined as below: 



We have modified the similarity transformations used by Liu 
et al. [9] for the current problem. Using (6), Eq.(l) is identically 
satisfied and Eqs. (2)-(5) become 



/"' + 



«+l 



if + g)f"-nif'+g')f' = Q, 



Pr 



m =g(0) = 0,/'(0) = 1, g'{Q)=A, 0(0) = 1, 



(7) 



(8) 



(9) 



(10) 



-y iv.= /"+>')g"\ , (11) 



r =^ r =^ 



where izx ^ind T;,, are the wall shear stresses and (y,,. is the wall heat 
flux. These are given by 



/ 8u dw\ ( 5v dw 



(12) 



using Eqs. (6) and (12) in Eq. (1 1), one obtains 



Rey2Q,.=/"(0), 



Rer'/27VM= - 



= /(0), 

-e'(0). 



(13) 



where Re^ = f7,v(A-+}')/v and Re^, = K,t,(A:+ j)/v are the local 
Reynolds numbers along the x — and y— directions respectively 
(see Table 1). The vertical component of velocity at the far field 
boundary can be expressed as 



w{x,y, 00 ) = — \Jav{x^yf ' 



n+1 



\f{<x>) + gi^)]. (14) 



/(«)^O,g(c»)^O,0(«)^O, 



Optimal homotopy analytic solutions 

The non-linear differential equations (7)-(9) with the boundary 
conditions (10) have been solved by optimal homotopy analysis 
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Figure 2. Comparison of analytical and numerical solutions for the temperature distribution. Lines: 15"^-order OIHAiVI solutions, Circles: 
Numerical solution. 

doi:1 0.1 371 /journal.pone.01 07287.g002 

method (OHAM) [34,35]. For this purpose, we first select tlie and the auxiliary linear operators are selected as below 

initial guesses /o, go and 9o of /()?), g(r]) and 6(t]) as under: 

^^■(^)=^-^' ^^(^)=rf;^-^' ^«(^)=^-^-(i^) 



/o('7) = l- exp(-i;), go('?) = ^[l - exp (-»;)], 
0o(i?)= exp (-)/), 



(15) 



If ^£[0, 1] is the embedding parameter and fi the non-zero 
auxiliary parameter, then the generalized homotopic equations 




Figure 3. Influence of stretching rates ratio ?. on the x— component of velocity /'. 

doi:1 0.1 371/journal.pone.01 07287.g003 
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corresponding to (7)-(10) can be written as follows 

{\-q)Lf[F{n- q)-M,l)]=qm,[F{rj- q),G{n; <?)], (17) 

(\-q)L,[G{n; q)-gM]=qhN,[F(n- q),G(ri; <?)], (18) 



N,[F(,r, q\ GOr, q)] 

d'G{,r,q) , 8'G(,r, q) 

(dF(n; q) dG(ir, q)\ SG{n; q) 



ill) 



8)1 J dr\ 



{\-q)Lg[0(n; q)-eom=qmg[FOr, ql G{n; q),0{n; q% (19) 



^('?;?)l,=o=o 

G(m q)\,,=Q=Q, 
SFin; q) 



SF(m q) 



drj 
8G(ir, q) 



(/=0 



(20) 



0, 



dt] 
dG{ r\; q) 



ri = 0 



dl1 



--A, 0(rj; g)L „ = 0, 



= 0, 0(,7; q) 



where the non-linear operators Nf, Ng and Ng are 
Nf[F{,r, q), G(n;q]\ 



8^F(ti; q) n+l 
8ri^ 1 



{F(n; q) + G(n; q)) 



s'Hm q) 

8ri^ 



(21) 



8F{)r, q) , 8G{n-, q)\ dF{ir, q) 



8r] 



dn J 



Ng[F{,r,q), Gim qmrnq)] 

1 8^0{ii; q) ( n+ 1 



Pr 8if ' V 2 

^d Fjir, q) ^ 8G(t]- q) 
drj 8r] 



{F(ii;p) + Gir,q)) 



q)- 



80(n\ q) 
dn (23) 



By Taylor's series expansion one obtains 

1 S"'F{,r,q) 



Hi; q)=Mn)+ E/^('7)r;/m07) = 



m\ dq'" 



, (24) 



q = 0 



Gin; q)=go(ti)+ '^gm(ti)q"';gm(ji)-- 



1 d"'G{tr, q) 



ml 8q"' 



00,; q)=e,{n)+ T e„,(>?)g"';e„,(»7)=^ '^"'?!'^^ 

^ ml dq'" 



q = 0 



q = 0 



, (25) 



(26) 




0 0.25 0.5 0.75 1 1.25 1.5 1.75 



Figure 4. Influence of stretching rates ratio ). on the y— component of velocity g'. 

doi:1 0.1 371/journal.pone.01 07287.g004 
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: -f"{0)____ 






n = 3 


I ^TMHjt^)] 








0.1 0.2 0.3 0.4 0.5 



0.7 0.8 0.9 1 1.1 1.2 



Figure 5. Influence of stretching rates ratio 1 on tlie si<in friction coefficients /"(O) and g"(0) and entrainment velocity /(x) + g(oo). 

doi:1 0.1 371/journal.pone.01 07287.g005 



Substituting ^ = 1 in the above equations yields the final result. 
The functions f,„, gm and 6„, can be determined fi-om the 
deformation of Eqs. (7)-(10). Explicitly the mth-order deformation 
equations corresponding to Eqs. (7)-(10) are as below 



i/[^,('?)-Z„/m-l('?)l=/<('/), 



(27) 



j3/- m-l 



k = 0 



n+l 

IT 



(31) 



drj 



dt] J dt] 



Lg[e,„{n)-X„A„-^{ri)]=hRi{il), 



(28) 



(29) 



m-l 



m-l 

dii 



n+l 



(fm-l-k +gm-l-k) 



k = 0 

m-l -A- , agni-\-k \ Clgk 



(32) 



dr\ J dr\ 



/m(0) = 0, 
gm(0) = 0, 



dfmin) 



nil) 



dt] 



i) = 0 



i) = 0 



=0, e,„(0)=o, 



(30) 



KM- 



1 c/^0m-l 

Pr dt]^ 



ni—l 

E 



k = 0 

m-l -A- , ^fem-l-A- 



«+l. .ddk 
2 l/m- 1 -A: +^m- 1 -k) 



(33) 



dii 



Ok 



sfmin) 

drj 



=0, e,„(oo)=o, 



Q,m<\, 
\,m>\. 



(34) 
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Figure 6. Influence of stretching rates ratio ). on the temperature 0. 

doi:1 0.1 371 /journal.pone.01 07287.g006 



In order to determine the optimal values of h we define the 
squared residuals of the governing Eqs. (7)-(10), C^, and as 



d - 











0 





/=0 y=0 



00 








0 





.,/=o /=0 



(35) 



M M M 

,/=o y=o ./=o 



(36) 



An. (37) 
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Such kind of error has been considered in other works [36-41]. 
The smaller Cm'-'' rnore accurate the mth order approximation 
of the solution. The optimal values of ft can be obtained by 
minimizing the Cji/'.?, through the command Minimize of the 
software MATHEMATICA (see Liao [36] for details). Alterna- 
tively MATHEMATICA package bvph 2.0 can also be used to 
calculate such values (see [41] for details). 

Numerical method 

Eqs. (7)-(9) subject to the boundary conditions (10) have been 
solved numerically by shooting method with fifth order Runge- 
Kutta integration procedure. First, we reduce the original ODEs 
into a system of order ODEs by substituting 

xi =f,X2 =f',X3 =f",x4=g,X5 =g',x6=g",xj = 6 and xs = 9' 
which gives 







Xj 








X4 




X5 




Xe 




x-i 




.-4 . 





X2 
X3 

(Xi +X4)X3+«(X2+X5)X2 
X5 
X6 

(Xi +X4)X6+«(X2+X5)X5 
Xg 

n+l 

Pr(^— (Xi +X4)X8-m(X2+X5)x7) 



n+l 



n+l 



Table 2. Numerical values of /"(O) and ^"(0) for different values of n and 1. 




n X 


/"{») 




g"(0) 






shooting 


bvp5c 


shooting 


bvp5c 


1 0 


-1 


-1 


0 


0 


0.5 


-1.224745 


-1.224742 


-0.612372 


-0.612371 


1 


-1.414214 


-1.414214 


-1.414214 


-1.414214 


3 0 


-1.624356 


-1.624356 


0 


0 


0.5 


- 1 .989422 


- 1 .989422 


-0.994711 


-0.994711 


1 


-2.297186 


-2.297182 


-2.297186 


-2.297182 



doi:l 0.1 371 /journal.pone.Ol 07287.t002 
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Table 3. Numerical values of local Nusselt number — 


6'(0) for various values of w, Pr and a. 




It Pr 












shooting 


bvpSc 


1 0.7 


0 


0.793668 


0.793668 




0.5 


0.972033 


0.972029 




1 


1.122406 


1.122321 


1 


0 


1 .000000 


0.999990 




0.5 


1 .224745 


1 .224742 




1 


1.414214 


1.414214 


7 


0 


3.072250 


3.072251 




0.5 


3.762723 


3.762724 




1 


4.344818 


4.344779 


3 0.7 


0 


1.292193 


1.292194 




0.5 


1.582607 


1.582607 




1 


1.827437 


1 .827427 


1 


0 


1 .624356 


1 .624356 




0.5 


1.989422 


1 .989422 




1 


2.297186 


2.297182 


7 


0 


4.968777 


4.968777 




0.5 


6.085484 


6.085485 




1 


7.026912 


7.026913 


doi:l 0.1 371 /journal.pone.Ol 07287.t003 



and the corresponding initial conditions are 



Xi ' 




- 0 " 


X2 




1 


Xi 




"1 


Xi, 




0 


Xs 






X6 




U2 


Xl 




1 


Xt, . 




.Ui. 



(39) 



Suitable values of the unknown initial conditions 
Ml =/"(0),M2 =^"(0) and K3 = 6'(0) are guessed and then integra- 
tion is carried out. The values of U\ ,U2 and i/3 are then iteratively 
computed through Newton's method such that the solutions satisfy 
the boundary conditions at infinity (given in Eq. (10)) with error 
less than 10"*. 



Results and Discussion 

This section contains the physical interpretations of the 
behavior of the interesting parameters entering into the problem. 
We compare the 15*-order OHAM solutions for temperature 6 
with the numerical ones for different values of «. Fig. 2 shows that 
data retrieved from both solution methods are identical, demon- 
strating the validation of our findings. 

Figs. 3 and 4 show the variations in horizontal and vertical 
components of velocity with an increase in stretching rates ratio I. 
It is clear that increase in k corresponds to an increase in the 
stretching rate along the >' — direction. Due to this reason the 
vertical component of velocity increases with an enhancement in X 



while the velocity in the x — direction decreases correspondingly. 
The wall velocity gradients f"(0), g"{0) and entrainment velocity 
f{co) + g{co) as functions of stretching rates ratio /. have been 
plotted in Fig. 5. Due to the bi-directional stretching sheet, there 
will be downward flow in the vertical direction. The vertical 
component at far field boundary is therefore expected to be 
negative in this situation. We notice that shear stresses at the wall 
increase when X is increased. Further, the larger values of X 
enhances the velocity of the cold fluid at the ambient. As a 
consequence, the entrainment velocity is an increasing function of 
X. 

Fig. 6 indicates that temperature 6 decreases with an increase in 
stretching rates ratio / for unity Prandtl number. Physically, an 
increase in X enhances the intensity of colder fluid at the ambient 
(as noticed in Fig. 6) towards the hot sheet which eventually 
corresponds to decrease the local fluid temperature. Fig. 7 
perceives the behavior of Prandd number Pr on the temperature. 
A bigger Prandtl number fluid has less thermal diffusivity and 
hence it allows less thermal effect to penetrate deeper into the 
fluid. As a result, temperature decreases and the thermal boundary 
layer becomes thinner when Pr is increased. This decrease in 
thickness of the thermal boundary layer is compensated with a 
larger wall slope of temperature function. 

Fig. 8 plots the wall temperature gradient against Pr with the 
variation in stretching rates ratio X. The wall heat transfer rate 
approaches the zero value for vanishing Prandtl number Pr — »0, a 
fact that is clear from the energy equation (9). Moreover, this Fig. 
compliments the results of Fig. 4. In bigger Prandti number fluids 
the convection is effective in transferring energy from the 
stretching sheet compared to pure conduction. Due to this reason 
the wall heat transfer rate is an increasing function of Pr. The 
reduction in thermal boundary layer thickness with an increase in 
X meets with the bigger magnitude of local Nusselt number. In 
other words the enhanced intensity of cold fluid at the ambient 
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towards the hot fluid closer to the sheet results in larger heat 
transfer rate at the sheet. 

Tables 2 and 3 provide the numerical values of skin friction 
coefiicients and local Nusselt number for different values of 
parameters by employing shooting method. The results are 
compared with the MATLAB built in function b\p5c and found 
in excellent agreement. We notice that wall shear stresses increase 
with an increase in X more rapidly at n = 3 when compared with 
n=\. The thinner thermal boundary layer accounted for larger n 
accompanies with larger temperature gradient along the sheet. 
The magnitude of increase in wall temperature gradient 0'(O) with 
an increase in Pr increases when n is increased. 

Conclusions 

For the first time, the flow and heat transfer over a plane surface 
stretched non-linearly in two lateral directions have been 
investigated. The simulation in this study assumes that the 
temperature across the sheet is non-linearly distributed. Both 
analytic and numerical solutions are obtained and found in 
excellent agreement. Following are the major results of this study. 
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